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Abstract 



A shadow of the moon, with a statistical significance of 5a, has been observed 
in the underground muon flux at a depth of 2090 mwe using the Soudan 2 de- 
tector. The angular resolution of the detector is well described by a Gaussian 
with a sigma < 0.3°. The position of the shadow confirms that the alignment 
of the detector is known to better than 0.15° and has remained stable during 
ten years of data taking. 

PACS 96.40.Cd, 13.85.Tp, 96.40.Tv 
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Typeset using REVT^X 



I. INTRODUCTION 



The Soudan 2 detector is a 963 tonne iron and drift-tube sampling calorimeter, designed 
to search for nucleon decay |3]||, situated in the Soudan mine in northern Minnesota, USA. 
The spatial resolution of the detector is < 1 cm in three dimensions. Throughgoing muons 
are typically reconstructed with track lengths of a few metres and the angular resolution for 
such tracks is expected to be a fraction of a degree. Because of the good angular resolution 
and the uniform exposure of the detector over the period of ten years from January 1989, 
the detector is excellently suited to a study of the possible existence of point sources of 
underground muons. 

A prerequisite for such a search is a thorough understanding of the alignment and angular 
resolution (point spread function or psf) of the detector. The angular diameters of the 
moon and sun are both ~ 0.5°. Both objects will occlude the high energy cosmic ray 
primaries responsible for the underground muon flux. Since their angular sizes are small 
and comparable with the psf of the detector the observation of shadows in the underground 
muon flux will confirm both the alignment and the angular resolution of the detector. 

The potential observation of moon and sun shadows with high energy cosmic rays was 
originally suggested by Clark ||; the observation of the effect of the solar magnetic field 
on the cosmic ray shadow of the sun has been suggested as a means of determining the 
composition of primary cosmic rays with energies of ~100 TeV Shadows of the moon 
and sun have been observed in large EAS arrays |§-§[. One experiment, sensitive to primaries 
with a median energy of 17 TeV, has observed a time-dependent displacement of the shadow 
of the sun which was attributed to the interplanetary magnetic field [PlJlOl . The observation 



of a moon shadow in the underground muon flux has recently been reported [ 11 , 12 1 . The 



data reported here were accumulated over a period of ten years during which the solar and 
interplanetary magnetic fields varied considerably and results are presented only for the 
shadow of the moon. The results of an analysis of the data for the shadow of the sun will 
be presented in a separate publication JT3[ . 

A shadow of the moon in the underground muon flux should be observable if the trans- 
verse momentum impulse, Ap t , imparted to a primary by the geomagnetic field is sufficiently 
small compared with the momentum of the primary. The Soudan 2 detector is located at a 
depth of 2090 mwe at latitude 47.82° N, longitude 92.24° W in NE Minnesota, USA. The 
local overburden is flat. The magnitude of the local geomagnetic field is 59 /xTesla and, 
because Soudan is close to the geomagnetic pole, the field at the detector is inclined at 
14.9° to the vertical, points down and to the north, and lies almost in the plane of the local 
meridian. An estimate, based on an impulse approximation and assuming a dipole field, 
of the transverse momentum imparted to a primary cosmic ray yields Ap t ~ cR^f3 A B 
where R® = 6378 km is the radius of the earth. When averaged over azimuth for primaries 
from directions close to the moon, Ap t is estimated to be ~ 25 GeV/c to the west. The 
p t impulses in the N — s and vertical directions average to zero. 10 TeV primaries would 
therefore suffer mean deflections of ~ 0.14° to the west, leading to their mean directions 
reconstructed in the detector being apparently displaced to the east. The rms value of Ap t 
is 50 GeV/c which would produce a broadening of the shadow of the moon of ~ 0.28° for 
10 TeV primaries. A more precise estimate of geomagnetic effects is described below. 

The position of the moon lies within approximately ±5° of the ecliptic which is itself 
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inclined at an angle of 23.5° to the equatorial plane. The declination of the moon therefore 
varies between ±28.5° throughout the year and its minimum zenith angle at Soudan is 
~ 20°. A muon arriving at the detector from a direction close to the moon must therefore 
penetrate at least 2200 mwe of rock to reach the detector, corresponding to a minimum 
surface energy of ~ 0.8 TeV. The mean slant depth is ~ 3000 mwe and the median surface 
energy of muons arriving underground from this direction is estimated to be ~ 1.5 TeV. 
Since the energy of a muon in an air shower is approximately one tenth of the energy of 
the primary, most muons observed in the detector originate from primaries with energies of 
approximately 6 to 15 TeV or greater when the moon is closer to the horizon. Geomagnetic 
effects are therefore expected to be small (O(0.1°)), if not entirely negligible, compared with 
the expected angular resolution of a fraction of a degree. 

Two other effects will influence the position and size of any shadow observed in the 
detector. Multiple Coulomb scattering of the muons in the rock overburden will broaden 
the psf of the detector and, since an object with fixed celestial coordinates moves in de- 
tector coordinates during a day, any misalignment of the detector will also contribute to a 
broadening of the psf as well as displacing the position of a shadow flT4 |. 

The broadening of any shadow depends upon the momentum spectrum of the primaries. 
The Monte Carlo simulation described in reference [|I5[] has been used to study the contri- 
butions to the broadening from geomagnetic deflections and multiple Coulomb scattering. 
A sample of 14,000 muons from cosmic ray primaries with initial directions within 3° of the 
position the moon (uniformly sampled throughout the period of data taking) was generated 
using the hemas cascade and sibyll |l7j hadronic interaction codes. The muons were 
tracked with geant to the Soudan depth. 

The mean energy of the primaries from these directions which resulted in a muon at the 
detector is 19 TeV. The impulse approximation described above, together with the energy of 
the primary, is used to estimate the geomagnetic deflection of the primary. The broadening 
due to multiple scattering (and the shower production mechanism) was determined by cal- 
culating the angle between the directions of the primary and the muon underground. The 
angular distributions (dN/d6 2 ) due to geomagnetic dispersion and multiple scattering are 
very non-Gaussian, being sharply peaked towards zero and having tails extending beyond 4 
degrees squared. The distributions are summarised in table | by the 50-percentile and mean 
values of 2 . The broadening due the geomagnetic field is small. Multiple scattering gives a 
mean square broadening of 0.43 degrees squared although the distribution is peaked much 
more closely towards zero than a Gaussian distribution with a = 0.46° which has the same 
< 6 2 >. The moon itself may be effectively taken as a point object. 

The angular size of the moon which will be deduced from the data will be characteristic of 
the psf of the detector for a point source where the psf folds together the angular dispersion 
due to detector resolution and alignment, geomagnetic dispersion and multiple Coulomb 
scattering. The detector angular resolution is expected to be a few tenths of a degree from 
considerations of track length and spatial resolution. 



II. DATA SAMPLE 

The Soudan 2 detector is described fully in references U and 0. It is composed of 
two layers of 1 x 1.1 x 2.5 m 3 modules containing corrugated steel plates interleaved with 
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1 m long, 15 mm diameter drift-tubes filled with argon/15% carbon dioxide at atmospheric 
pressure. The three coordinates of a track crossing a tube are obtained from an array of 
crossed anode wires and cathode pads situated at the ends of the tubes, and drift times. 
Pulse height, which is proportional to the ionisation deposited, is also recorded for each hit. 
The electron drift velocity is ~ 0.6 cm//xsec; the maximum electron drift distance is 50 cm. 
The detector has been surveyed to 0.1° about all three axes which are the vertical and very 
close (within 0.3°) to E - w and N - S. 

The detector has operated continuously since January 1989. Its size increased steadily 
to reach its present size of 14 m x 8 m x 5 m in late 1993. Data taking has been very stable 
with down time for maintenance and installation confined to daytimes. Muons are currently 
accumulated at a rate of approximately 7.5 million per year. A total sample of 58.5 million 
muons has been accumulated during the ten years until December 1998. 

Two complementary software algorithms are used to reconstruct tracks in the detector. 
The first, 'fmr', is a fast algorithm which assembles small, clean, track segments and ignores 
noisy regions due to catastrophic energy losses from bremsstrahlung or pair production. The 
second, 'search', associates hits in roads to form complete tracks. Both algorithms work 
with hits in the anode-time and cathode-time projections before a full three-dimensional 
track is defined. The electron drift velocities used by the two algorithms are independent 
and are corrected for changes in atmospheric pressure (which can amount to a few percent) 
and gas composition. In most cases a penetrating track will be reconstructed by both fmr 
and search. The search algorithm usually reconstructs a greater length of track but has 
somewhat worse systematics than fmr. The systematic difference between the parameters 
of tracks reconstructed by both search and fmr is < 0.05°. Both algorithms tend to be 
unreliable for tracks lying close to the principal planes of the detector. 

Stringent cuts are applied to the reconstructed data to ensure a high quality sample. 
The data are taken in runs which last approximately one hour. Any run which shows an 
anomalously high trigger rate or rate of reconstructed tracks, or other symptom of detector 
malfunction, is rejected. These cuts reject approximately 6% of the data runs. The track 
selection requires either a search track longer than 3 m which lies no closer than 10° to 
the principal planes of the detector or, if such a track is not found, an fmr track between 
1 m and 5 m long. In either case the selected track is also required to have a pulse height 
per unit length (i.e. dE/dX ) consistent with a minimum ionising particle and to have 
reconstructed local angles such that it would have triggered the detector. These cuts reduce 
the total sample of 58.5 million muons to a sample of 33.5 million for further analysis, with 
the major loss being the requirement on track length. 

Because of its good spatial resolution, and the many (typically ~ 100) hits used to form 
tracks, the inherent angular resolution of the detector is expected to be a small fraction of a 
degree for tracks a few metres long. The angular resolution can be investigated by examining 
the distribution of space angle, ip, between pairs of muons observed simultaneously (i.e. 
belonging to the same parent shower) in the detector. Figure [l] shows the distribution of 
if) 2 for a sample of 1.26 million muon pairs. For a Gaussian psf dN/dip 2 oc exp(— tp 2 /4a 2 ) 
where a is the angular resolution of the detector[|. The distribution shown in Figure p] can be 



4<r 2 appears in the denominator because the distribution for pairs of muons results from the 
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described by a Gaussian with a = 0.32° for ijj < 1°; the excess beyond 1° can be attributed to 
multiple scattering of the muons in the overburden. There is no evidence for any anisotropy 
in two dimensions. 



III. DATA ANALYSIS 
A. One dimensional analysis 



The time (utc) of each selected event is recorded with a clock synchronised to the wwvb 
time standard and used to convert the direction of the reconstructed track in detector 
coordinates to geocentric celestial coordinates (right ascension, a M , and declination, 5 M ). 
The celestial coordinates of the moon, a m and 5 m , at the event time are also calculated. 
The moon's position is corrected to topocentric coordinates using the latitude and longitude 
of the detector to account for the parallax due to the moon's relative proximity. The routines 
|18| used to calculate the lunar coordinates are accurate to better than 0.01° over the entire 



period of the exposure; the parallax corrections amount to as much as ±0.7° and —0.8° in 
right ascension and declination respectively. 

The space angle, 8, between the direction of the muon track and the direction of the 
moon is histogrammed. The angular density of muons, dN^/dO 2 , as a function of 9 is derived 
from this histogram by dividing the bin contents by the solid angle subtended by each bin. 
Figure || shows {l/n)dN^/ d6 2 versus 6. In the absence of a shadow this distribution should 
be flat. There is, however, a clear deficit of events in the first few bins which is attributed 
to the shadowing effect of the moon. The significance of the shadow is tested by comparing 
the difference in x 2 between the best fit to a flat distribution (x 2 = 82.9) and the best fit 
(x 2 = 58.3) to the form 

,1X " -nX(l-n(R 2 m /2a 2 )eM-0 2 /2a 2 )) (1) 



de 2 

where R m = 0.26° is the mean angular radius of the moon and A and a are free parameters 
representing the angular density of muons and the detector angular resolution respectively. 
The fitted function treats the moon as a point object at 9 — which removes XnR^ muons 
from the sample. Numerical studies indicate the the effect of the finite angular size of the 
moon is relatively unimportant if the detector angular resolution is of the order of 0.3°. 
The best fit values of the parameters are A = 607.1 ± 3.48 muons per square degree and 
a = 0.333° ± 0.048°. The value of A implies that a total of 128.9 muons (about one per 
month) are missing from the data sample due to the presence of the moon. The improvement 
in x 2 between the two fits of 24.6 for two degrees of freedom implies that chance that the 
observed deficit is a random fluctuation is < 10~ 5 : a shadow of the moon is observed with 
a statistical significance of 5er. 

Figure |2] clearly demonstrates the existence of a shadow of the moon in the underground 
muon flux and indicates that the psf of the detector is the expected small fraction of a 
degree. There is no evidence for any large misalignment of the detector. 



convolution of two similar distributions. 
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B. Two dimensional analysis 



The procedure to extend the analysis to two dimensions is straightforward and involves 
three stages. In the first, maps (two dimensional histograms) of the angular muon density 
in a region centred on the moon are constructed from the arrival directions of the muons. In 
the second, the maps are smoothed (rebinned) with a two-dimensional kernel to form a map 
of the local mean angular density of muons. The expected angular density in the absence 
of a shadow is estimated at the same time. Finally the expectation is subtracted from the 
smoothed data and Z y a test statistic which represents the normalised excess or deficit at a 
point, is evaluated. 

The celestial coordinates of a muon track are transformed into a coordinate system 
centred on the moon by first a rotation of — a m about the polar axis followed by a rotation 
of — 5 m around an axis in the equatorial plane. In this new system the x-axis points towards 
the moon and the y-axis lies in the equatorial plane. The x direction cosine, c x , is therefore 
equal to cos ip where ip is the angle between the track and the direction of the moon. To a first 
approximation c y ~ (a M — a m ) cos5 m and c z ~ 5^ — S m . For tracks at small angles (< 10°) 
to the direction of the moon, 1 — c 2 x = c 2 + c 2 ~ ip 2 , and the space angle, in degrees, is 
well approximated by ip = (180° /ti)(c 2 + c 2 )^ . Subsequent analysis is made using the quasi- 
angular coordinates 6*ew = (180°/7r)c 2/ and 6*ns = (180°/7r)c z which represent the angular 
distances between the direction of the muon and the moon in the west-east and south- 
north directions respectively. 

Tracks pointing towards the moon, i.e. with c x > 0, are selected. 2.44 x 10 5 tracks 
lie within an area of ±10° x ±10° around the moon. The coordinates #ew and #ns of the 
muons are used to make low- and high-resolution maps (two-dimensional distributions), 
covering 20° x 20° and 8° x 8° regions with the moon at their centres. The bin sizes used are 
A0 EW x A9 NS = 0.1° x 0.1° and A9 EW x A# NS = 0.04° x 0.04°. The 9 EW and # NS projections 
of the low resolution map are shown in Figure The flatness of the #ew distributions is 
a result of the many cycles of averaging in the east-west (right ascension) direction during 
the long exposure of the detector. The rise in the number of muons with #ns (declination) 
results from the decreasing overburden for muons with small zenith angles. (The moon is 
never overhead at the detector latitude of 47.8°.) 

The density, A(# EW , 9 NS ), of muons in the maps expected in the absence of a shadow 
is estimated by assuming that A factorises, in other words that the number of events 
expected in a bin centred at #ew, #ns can be represented by dN^ = Ag^ew^ns = 
F(^ew) x ^(^ns)^ew^ns- The functions Y(9ew) and Z(6*ns) are obtained by fitting 
quadratics to the 9 EW and 9 NS projections. The results of the fits for the low resolution 
map are shown in Figure |3|. The fits are good with % 2 /ndf = 189.1/198 for the #ns projec- 
tion and x 2 /ndf = 170/199 for the 6*ew distribution. 

The fitted functions Y(9ew) and Z(#ns) are used to predict the muon density, 
A(^ew,6 i ns) = d 2 N ll /d9Ewd9^s, in the absence of a shadow at any point in the map. As 
can be seen from Figure |3], A varies by a factor of three over the range of —10° to +10° in 
6*ns direction but can be taken to be flat on a scale of 1°, allowing quantitative predictions 
to be made for the two-dimensional analysis which is described below. For the both the 
low- and high-resolution maps, the x 2 between the number of events in any bin, riy, and the 
prediction, my = Y(9 EW ) x Z(# ns )A# E wA#ns, is satisfactory with rV nd f = 40076.6/39994 
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for the coarse map and xV n df = 40559.4/39994 for the fine map. 

The predicted angular density of muons at the centre of the low resolution map is 
A(0,0) = 608.3 muons per degree squared, consistent with the angular density obtained 
from the fit to the one-dimensional dN^/dO 2 distribution. The method used to estimate A, 
however, includes the central region where there is a known deficit of events (caused by the 
shadow of the moon) and therefore results in a slight underestimate of the expectation. The 
underestimate will be greater for the finer binned map; the predicted angular density of 
muons at the centre of the high resolution map is A(0, 0) = 592.1 muons per degree squared. 
To avoid any bias caused by this underestimate, the value of A is taken uniquely from the 
low resolution map in the subsequent stages of analysis. 

A two-dimensional Gaussian kernel of angular width is used to rebin the data in the 
maps. The contents, nij, of bin i,j are replaced with dij = J2kJ2i n (i+k)(j+i)Wki- The weight 
Wki = (l/2cr^) exp(— ipl l /2al)A9 EW A9 NS , where ipki is the angular distance between bin i, j 
and bin k, I. The sums over k and / include all bins within 3<Tfc of bin i, j and J2k J2i w m — 1- 

The expected number of events, rriy, in bin i,j in the absence of a shadow, is predicted 
from the fits to the projections of the map described above. The resultant map of rriij 
is then treated identically to the raw data map to obtain the rebinned expected angular 
density b^. Maps are then made of the normalised statistic Z^ = (d^ — b^)/ y / var(6jj) where 
var(&jj) = J2k J2i m (i+k)(j+i)wli- The value of A is sufficiently large that this procedure, which 
treats the statistical fluctuations of the background as Gaussian, introduces negligible error 
despite the inherent Poisson fluctuations of the small number of events (typically eight) in 
any bin [TjJ . 



Since the bin size of the maps is small compared with a (the angular resolution of the 
detector) and (the rebinning kernel), and the background is flat over a region ^> cr, the 
expected value of Z for a point source (or sink) of strength N m may be estimated by treating 
the summations as integrals to yield 



4 \ N m (2) 



and 



[al + a 2 ) a k 



(3) 



Z is maximised when = cr, i.e. when the rebinning kernel is matched to the angular 
resolution of the detector. Equation (|j) predicts Z m = —4.43 for the shadow of the moon 
using A(0,0) = 608.3 muons per degree squared and matching ak to the value cr = 0.333° 
found from the fit to the dN^/ d9 2 distribution described in section |111 A| . 

A contour map of normalised deviations, Z, made with a kernel size of = 0.333° 
for a ±7° x ±7° region centred on the moon is shown in Figure |j. To avoid confusing 
detail, contours are shown only for the regions where \Z\ > 2.0. The shadow of the moon 
is unambiguously visible as the deep minimum [Z = —4.69) at the centre. The 25 extrema 
with \Z\ > 2.0 in the entire area of 196 square degrees are consistent with the 31 expected 
from statistical fluctuations 1 12011 ; there are no other minima with Z < —3.0. The observed 



value of Z = —4.69 for the shadow suggests that the angular resolution is somewhat better 
than 0.333° determined from the one-dimensional analysis. 
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The bin size of 0.1° x 0.1° is rather too coarse compared with the angular resolution 
of approximately 0.3° to allow a precise determination of the position of the shadow to be 
made. A further analysis is made using the finer binned map with a bin size of 0.04° x 0.04°. 
Figure [| shows a contour map of Z for a ±2.8° x ±2.8° region using a kernel size of a^ = 0.29° 
which was chosen to minimise Z for the shadow. In this map A has been evaluated using 
the parameters determined from the fits for Y(9ew) and Z(#ns) to the coarser map to avoid 
the underestimate discussed above. 

With the finer binning, the minimum, which has a depth Z = —4.98, at (#ew = 
0.1°, #ns = 0.1°) is within 0.15° of the nominal position of the moon. The Z = —4.5 
contour, which corresponds to ~ 75 % CL, encompasses the origin and there is therefore 
no evidence that the detector is incorrectly aligned. It is not possible to separate what 
fraction of the small observed displacement of the shadow could be due to misalignment, 
the statistics of the sample or geomagnetic effects, although the displacement of 0.1° east is 
consistent with the estimated geomagnetic displacement of 0.076° discussed in section |. 

IV. LONG TERM STABILITY OF THE SOUDAN 2 DETECTOR 

As described in section [I] the Soudan 2 detector was assembled over a period of five years 
ending in 1993. Some modules were removed from the detector and refurbished in 1994 and 
1995. The detector continued to operate during these periods and was re-surveyed after 
each movement. The clean shadow of the moon demonstrated in sections |111 A| and |111 B| , 
and the statistics of the sample allow the data to be subdivided to check the stability of the 
alignment of the detector. 

The data have been divided into two approximately equal samples for the periods of 
1989 to 1994 inclusive, and 1995 to 1998. The analysis described in section |Tl rB| has been 
repeated for each period. The high resolution maps, made with au = 0.29°, for both periods 
are shown in Figure || A shadow of the moon is clearly seen in each map. The positions 
and depths of the shadows are given in table The depths of the minima are both close 
to the expected value of Z = — 4.98/a/2 = —3.52 expected for half the total sample. The 
positions of these shadows confirm that the alignment of the detector has remained stable 
during the ten year period. 

V. CONCLUSIONS 

An shadow of the moon in the underground muon flux has been observed using the 
Soudan 2 detector with a statistical significance of 5a. The size and position of the shadow 
confirm that the alignment of the detector is correct to better than 0.15° and that the psf of 
the detector can be adequately described by a Gaussian with a = 0.29°. The geomagnetic 
deflection of the primaries, which is expected to be of the order of 0.1° prohibits a more 
precise statement about the accuracy of the alignment of the detector. The subdivision of the 
data into two periods confirms that the alignment of the detector has remained stable over 
the ten years from 1989 to 1998. These data confirm the excellent capabilities of Soudan 2 
as a detector of possible astrophysical point sources of underground muons. 
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TABLES 



Source of broadening # 2 < 6 2 > < 8 > 

(degrees 2 ) (degrees 2 ) (degrees) 

Moon (flat disc) 0.034 0.034 

Geomagnetic 0.04 0.32 0.076 west 

Multiple scattering 0.095 0.43 0.0 

Gaussian (cr = 0.46°) 0.3 0.43 

TABLE I. Calculated contributions to the broadening of a shadow of the moon. The first 
column gives the 50-percentile of the distribution; the last row gives the values for a Gaussian 
distribution for comparison. The last column gives the mean displacement of a particle trajectory; 
the reconstructed trajectory, and hence any shadow, is displaced in the opposite sense. 



Period A(0,0) Z min 8 E w #ns 

1989 - 1998 608.3 -4.98 0.1° 0.10° 

1989 - 1994 295.9 -3.78 0.1° 0.02° 

1995 - 1998 311.8 -3.39 0.1° 0.18° 

TABLE II. Observed positions of the shadow of the moon for the entire data sample, and the 
periods 1989 - 1994 and 1995 - 1998. 
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FIGURES 
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FIG. 1. The distribution, dN/dijj 2 , of the square of the space angle between pairs of muons 
observed simultaneously in the detector. The line is the result of a fit dN/dip 2 oc exp(— ip 2 /4a 2 ) in 
the region -ip 2 < 1 degree 2 . The fitted value of a is 0.32°. 
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FIG. 2. The angular density of muons, (l/^dN^/dO 2 , versus 9, the angular distance between 
the muon direction and the nominal position of the centre of the moon. 
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FIG. 3. The 6ew and 6*ns projections of the raw data in the 20° x 20° low resolution map. 
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FIG. 4. Contour map of normalised deviations, for a ±7° x ±7° region centred on the moon 
made with a rebinning kernel with ct/% = 0.333°. The contour lines are spaced by AZ = 0.5 and 
shown only where \Z\ > 2.0. Regions with Z > are shaded. 
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FIG. 5. Contour map of normalised deviations, Z, for a ±2.8° x ±2.8° region centred on the 
moon made with a rebinning kernel with = 0.29°. The contour lines are spaced by AZ = 0.5 
and shown only where \Z\ > 2.0. Regions with Z > are shaded. 
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FIG. 6. Contour maps of normalised deviations, Z, for a ±2.8° x ±2.8° region centred on the 
moon made with a k = 0.29° for the periods 1989 - 1994 (left) and 1995 - 1998 (right). The contour 
lines are spaced by AZ = 0.5 and shown only where \Z\ > 2.0. Regions with Z > are shaded. 
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